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ABSTRACT 

When an open system of classical point particles interacting by Newtonian gravity 
collapses and relaxes violently, an arbitrary amount of energy may in principle be 
carried away by particles which escape to infinity. We investigate here, using numer- 
ical simulations, how this released energy and other related quantities (notably the 
binding energy and size of the virialized structure) depends on the initial conditions, 
for the one parameter family of starting configurations given by randomly distribut- 
ing A'' cold particles in a spherical volume. Previous studies have established that the 
minimal size reached by the system scales approximately as N~^^^ , a behaviour which 
follows trivially when the growth of perturbations (which regularize the singularity of 
the cold collapse in the — > oo limit) are assumed to be unaffected by the bound- 
aries. Our study shows that the energy ejected grows approximately in proportion to 
TV^/'^, while the fraction of the initial mass ejected grows only very slowly with N, 
approximately logarithmically, in the range of N simulated. We examine in detail the 
mechanism of this mass and energy ejection, showing explicitly that it arises from the 
interplay of the growth of perturbations with the finite size of the system. A net lag of 
particles compared to their uniform spherical collapse trajectories develops first at the 
boundaries and then propagates into the volume during the collapse. Particles in the 
outer shells are then ejected as they scatter through the time dependent potential of 
an already re-expanding central core. Using modified initial configurations we explore 
the importance of fluctuations at different scales, and discreteness (i.e. non-Vlasov) 
effects in the dynamics. 

Key virords: Virialization; spherical collapse; TV-body simulations 



1 INTRODUCTION 

That isolated systems of initially cold self-gravitating parti- 
cles collapse and relax violently to produce virialized quasi- 
equilibrium structures has been known for many decades — 
essentiall y since numer ical simulation of such systems began 
(see, e.e:.. lHenonl (| 1964 )1. In cosmology it has emerged in the 
last decade or so, through large numerical simulations, that 
a good approximate description of the structures formed 
can be given in terms of "halos", which are similar quasi- 
equilibrium structures formed from collapse of matter in fi- 
nite regions, containin g initially cold ( dark matter) particles 
(for a review, see e.g. ICoorav fc Sheth (2002)). Analytical 
understanding of the strongly non-linear physics involved in 
the formation of these quasi-equilibrium states is, despite 
the importance of the problem and many attempts to solve 
it, extremely limited. For example, while the simple density 



profiles which result from a class of uniform initial condi- 
tion s have been known since early numerical studies (see, 
e.g.. Ivan Albadal l|l982l) and references therein), there is still 
no theory which can clearly explain them. Likewise, in cos- 
mology, numerical evidence for the "universality" of certain 
simple halos profiles (difi'erent to those obtain ed from quasi- 
uniform col d collapse) has been produced (|Navarro et al.l 
(|l996l . ll997l )). but understanding of the physical reason for 
the ubiquity of these profiles is still lacking. 

In attempts to gain greater insight into the physics lead- 
ing to such virialized states, with the goal notably of under- 
standing the dependence of their properties on initial condi- 
tions, one angle of approach, which is the one adopted here, 
is to study in detail the evolution from some limited class 
of initial conditions. We focus here on what appears to us a 
neglected, and potentially important, aspect in the study of 
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collapse and virialization of open self-gravitating systems: in 
the phase of (violent) relaxation such systems may, in gen- 
eral, eject a fraction of their mass which carry away a 
finite amount of energy as kinetic energy, i.e., some of 
the particles come out of the collapse with positive energy 
and may escape to infinity. The binding energy and charac- 
teristic size of the virialized state, in particular, are directly 
related to these ejected quantities. Indeed, using (1) total 
energy and mass conservation, (2) the virial condition for 
the bound particles, and (3) neglecting the potential energy 
between the bound and escaping particles, and that of the 
escaping particles themselves, it follows that 



W" = 2{Eo - K") 



(1) 



where W" is the (negative) potential energy of the bound 
particles, and Eo is the total initial energy. While for iSo > 
some energy must be ejected, for Eq < energy may, or may 
not, be ejected. In principle there is no upper bound on this 
ejected energy (as there is no lower bound on the gravita- 
tional potential). Indeed it is known that a self-gravitating 
system of point particles is intrinsically unstable towards 
the ejection of an infinite amount of energy — the so- 
called "gravo-thermal catastrop he" l|Lvnden-Bell fc WoodI 
Il968l : iBinnev fc Trema"in3ll994l ). This instability is under- 
stood, however, to be relevant in practice only on time scales 
which diverge with particle number A'^. The term "violent 
relaxation", on the other hand, refers precisely to a phys- 
ical relaxation process which takes places on much shorter 
scales, in volve purely mean field (i.e. "coUisionless" ) physics 
(|Lvnden-B cll 1967) . In the study we present in this paper we 
will see that the energy ejected on these shorter time scales, 
during violent relaxation, appears to be unbounded above 
for the class of initial conditions we study. Further we will 
investigate whether the ejection we observe in our simula- 
tions can be fully understood as a mean field phenomenon. 

The central considerations in this article — energy ejec- 
tion and the validity of the mean field limit in cold collapse 
— are of relevance in theoretical attempts to understand vi- 
ole nt relaxation . For e xample, the original theory proposed 
bv lLvnden-Beil (|l967l ) to predict the properties of virialized 
states made both the assumption that these states have the 
same energy and mass as the initial conditions, and that the 
dynamics is purely mean field (i.e. governed by the Vlasov- 
Poisson equations). Although this theory is believed to be 
inadequatqil — and in any case it may not be applicable to 
the particular set of cold initial conditions we study — such 
assumptions are usually made in theoretical approaches to 
understanding these states. It is thus important to under- 
stand the extent to which they are valid. Interestingly we 
find in our simulations that the shape of the density pro- 
files of the relaxed systems appear to be very robust, i.e., 
independent of the initial conditions, despite the fact that 
their global macroscopic parameters, notably their mass and 
energy, vary with A''. 

In this article we address these questions primarily us- 
ing the very restricted set of initial conditions given by dis- 
tributing A'^ particles randomly in a spherical volume, and 
ascribing zero velocity, i.e., a sphere of cold matter with 



Poissonian density fiuctuations. The simplicity of these ini- 
tial conditions is that they are characterized by the sin- 
gle parameter A''. On the other hand, they present an in- 
trinsic numerical complexity: while the problem is well de- 
fined for any finite A'^, it tends formally, as A'^ ^ oo, to 
the exactly uniform spherical collapse which leads to a 
density singularity at a finite time. This makes it expen- 
sive to integrate numerically for increasing A'^. Indeed for 
this reason many authors have excluded it from numeri- 
cal studies, focusing instead on spherical models with non- 
trivial inhomogeneous distributions (e.g. radially dependent 
densi t y) and /or significant non-zero vel o cities ly an Albad 
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McGlvnnI 1 19841: IV illumse nl 1 19841: lAguilax fc Merrit, 
Theis fc Spurzeml'll999i: .Merrall fc HenriksenI |2003| : 



Rov fc Pere2jl2004l : lBoilv fc Athanassoulall2006l l. As we will 

see, it is precisely the fact that the system collapses by a 
large factor, allowing particles to reach very large velocities, 
that leads to the very significant energy ejection which we 
focus on. We will discuss, in the part of the paper on the 
mean-field approximation, some other very specific initial 
conditions. We will also consider in our conclusions the ex- 
tent to which our study may be pertinent to other initial 
conditions. 

We note also that these initial conditions are the most 
evident discretisation of the exactly uniform spherical col- 
lapse model, which is a reference point in the theory of non- 
linear evolution of self-g ravitating system s in a strophysics 
and cosmology (see, e.g. ICoorav fc ShethI (|2002l )). The as- 
sumption usually made in using this model to predict the 
behaviour of physically relevant systems (e.g. in cosmology, 
via the formalism of Press and Schecter) is precisely that 
the total initial mass and energy are virialized. Our study 
shows that this is not a good approximation for this class of 
"almost" uniform spherical collapses. 

Despite the numerical difficulties associated with these 
initial conditions, several extensive studies of them have, 
however, been reported in the literature, most notably a de- 
tailed stud v bvlAarseth et al. (|l988l ). and a more recent one 
reported in lBoilv et al. (2002). We will discuss these works 
in greater detail below, and use some of their results. Both 
works focus on how the singular collapse of the infinite A'^ 
limit is regulated when there are a finite number of par- 
ticles, and in particular the scaling observed i n numerical 
simulations for the minimal size of the systerrQ. iBoilv et al.l 
(2002) also considers a wider class of non-spherical collapses, 
finding results to which we will return briefiy in our con- 
clusions. Our study can be seen as an extension of, and is 
complementary to, these studies, with which our results are 
in agreement. In more recent work llguchi et aL l (|2005l . l2006D 
have also studied collapse from such initial conditions, fo- 
cussing on the properties of the virialized state. Although 
they do not consider the ejected energy quantitatively, these 
authors do remark on its potential importance in open 
systems. Some of t he works ci t ed ab o ve on coUisionless 
relaxation (notably Ivan Albadal l|l982D : lAguilar fc Merritd 



^ For recent stud ies see e.g. lArad fc Johans sonI 

lArad fc Lvnden-Belll ||2005| '). and also that of lLevin et al.l 
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^ In the context motivating the study of lAarseth et al] ll 19881 ) 
the "points" are actually masses with extension (e.g. proto-stars) 
and the central question the authors wish to address is whether 
these masses survive or not the collapse of a cloud of which they 
are the constituents. 
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(| 19901 ') : IRov fc Perej l|2004 )') do include some quite cold uni- 
form spherical initial conditions, with initial virial ratios of 
order 0.1, and note that there is significant mass ejection 
for this case. For the relatively small number of particles 
considered in these cases the energy ejection is, however, 
modest and does not strongly modify the final state, and 
the dependence of this quantity on particl e number has not 
been explorejfj We note also the studies in lDavid fc TheunsI 
l|l989l ):lTh euns fc David ( 1990i) , which explicitly discuss and 
model the ejection of particles ("escapers") from pulsating 
spherical systems. 

We will describe here not only the dependence of the 
energies, and various other quantities, on A'^, but also detail 
the physical mechanism which leads to the ejection of mass. 
The probability of ejection turns out to be closely correlated 
with particles' initial radial positions, with essentially parti- 
cles initially in the outer shells being ejected. The reason for 
this correlation is simply that these particles which are ini- 
tially near the outer boundary systematically lag (in space 
and time) with respect to their uniform spherical collapse 
trajectories more than those closer to the central. These 
particles then gain the energy leading to their ejection in 
a very short time around the collapse time, as they pass 
through the time-dependent potential of the particles ini- 
tially closer to the centre, which have already collapsed and 
"turned around" . This lag, which has a very non-trivial de- 
pendence on the initial radial position, is manifestly a phys- 
ical effect coming from the boundary; fiuctuations to uni- 
formity evolve differently depending on their position with 
respect to the boundary. Indeed the lag which leads to the 
mass ejection first develops at the outside of the system and 
propagates into the volume during the collapse phase. The 
importance of the finite size of the system in this respect is in 
contrast to t he physics required i o und erstand the s caling s 
observed by lAarseth et al.l l|l988h and iBoilv et all (|2002l i. 
These can be understood by analyzing only the growth of 
fiuctuations in the collapsing system in the approximation 
that it is of arbitrarily large size, which corresponds to the 
case of a matter dominated contracting universe. Neverthe- 
less, as we will explain, we can use these arguments also to 
understand the scaling we observe of the ejected energy once 
the lag is assumed to be produced at collapse by the mech- 
anism desribed above. While this is highly consistent with a 
description of the ejection phase which is clearly mean-field 
like — the lagging particles propapating through the time 
dependent potential of the rest of the mass — it is not, as we 
discuss, evident whether this is true of the physics involved 
in the development of the lag (and therefore determination 
of the ejected mass). In the final section we address more 
precisely what this mean field limit is, and how one should 
extrapolate numerically towards it. We then report some nu- 
merical tests on appropriate modifications of the initial con- 
ditions, which probe such convergence. While we establish 
that there is indeed good evidence for convergence, there 



are still measurable fiuctuations due to non-mean field ef- 
fects in macroscopic quantities such as the ejected energy at 
the largest particle numbers we have simulated. 

The article is organized as follows. In the next section 
we recall the predictions for the A'' dependences of vari- 
ous quantities which follow from the uniform and perturbed 
spherical collapse model. In the following section we describe 
our simulations and principle results concerning the A'' de- 
pendence of the final state produced by violent relaxation. In 
Sect.|4]we investigate in detail the mechanism which leads to 
the ejection of mass and energy which we have observed. In 
Sect. Owe discuss the validity of the Vlasov-Poisson limit in 
describing the evolution we have studied. Finally we give a 
summary of our findings and conclusions in the last section. 



2 SCALINGS IN THE PERTURBED 
SPHERICAL COLLAPSE MODEL 

In this section we recall basic results on the exactly uniform 
spherical collapse model, and then the scalings of physically 
relevant quantities (minimum attained radius, time of col- 
lapse, velocities at collapse etc.) which are obtained from a 
straightforward analysis of the evolution of perturbations to 
this model. 



2.1 Uniform spherical collapse 

The radial position r{t) of a test particle in an (idealized) ex- 
actly uniform spherical distribution of purely self-gravitating 
matter of initial density po and initially at rest (at time 
f = 0) is simply given by the homologous rescaling 



r{t) = R{t)r{Q) 



(2) 



where the scale factor R{t) may be written in the standard 
parametric form 



R{i) = -{l + cos{0) 



(3) 



(C + sin(C)) 



and 



32Gpo ■ ^"^^ 

At the time Tscm the system collapses to a singularity. It will 
be useful in what follows to recall how the physical quanti- 
ties diverge close to this singularity. Taking ^ = vr — e and 



expanding to leading order in t gives {t — Tscm) 
which it follows that 



R(t) ~ - Ts. 



,2/3 



e , from 



(5) 



and therefore the test particle velocities v{t), proportional 
also to the initial radius r(0), scale as 



v(t) ~ [t - 



'1/3 



(6) 



^ lAguilar fc MerrittI ill 9901) actually discount a possible system- 
atic dependence on Af in a footnote on page 36. 
■* We note that this formation of a "core-halo" st ructure during 
collapse is remarked on and briefly discussed by lAarseth et al.l 
lll988l ). The link to mass/ energy eject i on are not, however, dis- 
cussed in this article, or in iBoilv et al.l ll2002h . 



2.2 Perturbed spherical collapse 

In this exactly uniform limit the evolution is independent 
of the size of the system, and the "particles" do not see its 
finite size (until the collapse). These solutions are thus for- 
mally valid as the radius of the system tends to infinity, and 
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indeed they are precisely those for the evolution of an infi- 
nite contracting universe containing only matter, which are 
known to coincide with those derived in general relativity. 
The last equations given above then correspond to the be- 
haviours in a contracting Einstein de Sitter (EdS) universe 
obtained when the curvature (corresponding to the initial 
cold start) is neglected. 

The system we study here — randomly placed par- 
ticles in a spherical volume — can be treated, up to some 
time and at sufficiently large scales, as a perturbed version 
of this uniform limit. While in general the associated per- 
turbations about uniformity will be expected to evolve in a 
way which will be sensitive to the finite size of the system, 
an approximation one can make is that such effects are neg- 
ligible, i.e. to treat the perturbations as if they evolve also 
in an infinite contracting system. Quite simply this means 
we neglect the effect of the boundaries on the evolution of 
the density perturbations. 

In the manner standard in cosmology (for the case of 
an expanding universe) one can then consider the fluid limit 
for the system a nd solve the appropriate equations pertur- 
batively (see e.g. IPeeblesI (jToSQ,)). In the eulerian formalism 
this gives, at linear order, a simple equation for 5(x), the 
density fluctuation (with respect to the mean density): 



S + 2HS - A-kGpoS = 



(7) 



where H(t) — R/R (dots denotes derivatives with respect to 
time) is the contraction ("Hubble") rate. These equations 
are derived in "comoving" coordinates x — v/R(t), where r 
are the physical vector positions. Note that 



i?x = r — i?x , 



(8) 



i.e., R{t)x{t) is the velocity of the particle minus the velocity 
it would have in the limit of uniformity (in the cosmological 
context, the "peculiar" velocity with respect to the "Hubble 
flow"). 

It is straightforward (see Appendix A) to solve Eq. (Q 
by rewriting is as an equation for 5{R). In the limit R <^ 1, 



S{R) ~ R 



Sit) ~ [t - 



-3/2 



(9) 



(10) 



This is simply the usual decaying mode of the expanding 
EdS universe, which becomes the dominating growing mode 
in the contracting case. 

A much more detailed analytic treatment of perturba- 
tions to the spherical collapse m odel in the fluid limi t, for 
the finite system, may be found in lAarseth" et al.l(| 19881 ). We 
give only the above results here, for the infinite radius limit, 
as they are the only ones we will make use of below. 



2.3 Predictions for scalings 

The singular behaviour of the spherical collapse is regulated 
by the fluctuations present at any flnite in the initial 
conditions we study. A simple estimate of the scale factor 
Rmin at which one expects the spherical collapse model to 
break down completely may be obtained by assuming that 
this will occur when fluctuations at some scale (e.g. of order 
the size of the system) go non-linear. For Poisson distributed 



partic les we have a mass variance (see e.g. iGabrielli et al.l 

Jiooi)) 



cr (r) 



{{N{r) ^ {N{r))r 



(11) 



{N(r)y (Nir)) ' 

where the angle brackets denote an ensemble average over 
realizations. Using this as the initial normalisation of the 
density fluctuations, and the growth given in Eq. ((Ojl, we 
can infer 



Rn 



oc iV" 



-1/3 



(12) 



Using the scalings in Eqs. (O and © it then follows that 
the maximum time tmax until which the spherical collapse 
model is a reasonable approximation is expected to scale as 



[imax 



OC N' 



-1/2 



(13) 



and the maximal infall velocity at any given radius as 

Vmax OC N^^^" . (14) 

Note that we do not need to make explicit any scale in 
the argument to obtain this result: the scaling with TV is ob- 
tained simply by taking the criterion that some amplitude 
is reached by the fluctuations at some fraction of the system 
size. In fact the result can be obtained on purely dimensional 
grounds: given that we are neglecting the finite size of the 
system, the only length scale in the problem — as gravity 
itself furnishes no scale — is given by the mean interparticle 
distance £ oc A^"^''"^. Clearly then the length scale deter- 
mined when Rmin is combined with the size of the system 
must scale in this way, if the physical processes determining 
it are ind eed independen t of th e size of the system. 

Both lAarseth et all (|l988l ) and iBoilv et al.l l|2002l ) use 
simple, but equivalent, variants of the above argument to 
obtain these same scalings. In particular they can be formu- 
lated in terms of the spread in the collapse time of approxi- 
mately spherical overdense and underdense regions, with ini- 
tial amplitude fixed by Eq. (|lip . Or, alternatively, one can 
analyse the scaling of the pressure associated with growth 
of the peculiar velocities which must compensate the inward 
gravitational pressur e to stop the c o llapse . One of the cen- 
tral fi ndings of both lAarseth et al.l l| 19881 ) and iBoilv et al.l 
(l2002h is that the scaling Eq. (|12p is in fact observed, and 
we will verify again in our simulations that this is indeed 
the case. 

The fact that these scalings are observed — and thus 
that the underlying approximation of (i) neglecting the 
boundaries of the system, and (ii) treating the fiuctuations in 
the fiuid limit, appears to work well — makes it instructive 
to compare collapsing spheres with different numbers of ini- 
tially Poisson distributed particles at appropriately rescaled 
times. Let us consider two such spheres, with say A''i and 
A^2 > A''i particles, as if they were just two spheres of dif- 
ferent radii (_Ri and R2 = Ri{N2/Ni)^^^) evolving in an 
infinite collapsing universe with initially Poisson distributed 
points. In the approximation that the perturbations grow 
as given by linear fiuid theory [i.e. as in Eq. ((UJ above], 
the sphere with N2 points should be exactly equivalent, up 
to transients due to the cold start, to that with A''i points 
at the initial time after the evolution of the scale factor 
brings its initially smaller density fiuctuation (of amplitude 
OC N2^^'^) to the initial larger amplitude (oc N^^^"^) of the 
sphere with TVi particles. This gives precisely the rescaling 
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inferred above, i.e., the different spheres should more gener- 
ally — and not just at the time corresponding to Rmin — 
be equivalent at scale factors rescaled in proportion to N^^^ . 
Deviations from equivalence in these rescaled time variables, 
on the other hand, must arise from physical effects not de- 
scribed by this approximation of neglecting the finite size 
of the system in its evolution up to collapse. We will thus 
compare below the different A'^ systems at these different 
scales factors, and discuss to what extent the observed scal- 
ings of the ejected mass and energy can be understood in 
this simple approximation. 



COLD SPHERICAL COLLAPSE: 
RESULTS FOR POISSONIAN 
FLUCTUATIONS 



BASIC 



We describe in this section numerical results for evolution 
from initial conditions in which TV particles are distributed 
randomly in a sphere, and given zero initial velocity. Be- 
sides reproducing known results for various quantities — 
notably the minimal collapse radius (or maximal potential 
energy), and the final profiles of the virialized structures 
— we present also our results for the TV dependence of the 
ejected mass and energy. 



3.1 Initial conditions and choice of units 

Table [1] shows the names of the different simulations and 
the associated particle number TV. The i-th realisation of 
an TV particle configuration is thus labelled PN — i. The 
parameter I is the mean interparticle separation in the initial 
configuration, defined as ^ = (SV/AivNy^^ = R/N^'^ where 
V is the volume of the sphere of radius R. The unit of length, 
here and throughout the paper, is the diameter of the initial 
sphere, i.e., R = 0.5. The parameter m is the (identical) mass 
of the particles. The unit of mass, here and throughout the 
paper, is chosen so that the initial mass density po is unity, 
i.e., Po = mN/V — QmN/n — 1. Finally we take our unit 
of time equal to Tscm, the uniform spherical collapse time as 
defined in Eq. @ Q. 

The parameter e is the softening parameter introduced 
in the numerical integration, which we will discuss at length 
below. 



3.2 Numerical code and parameters 

We have performed numerical simulations us- 
ing the pubhcly available code G AD GET2 
dwww ■ mpa-gaxching . mpg . de/ gadget /right . htmlj |2000| : 
ISpringellbOOsi r This code, which is based on a tree algo- 



rithm for the calculation of the gravitational force, allows 
one, in particular, as desired here, to perform simulations 
of a finite system with open boundary conditions. The 
two-body potential used is exactly the Newtonian potential 
for separations greater than the softening length e, and 



Name 


N 


£ 




e/e 


m 




P512-1 


512 


0.63E- 


01 


0.05 


O.lOE 


02 


P512-2 


512 


0.63E- 


01 


0.05 


O.lOE 


02 


P512-3 


512 


0.63E- 


01 


0.05 


O.lOE 


02 


PI 024-1 


1024 


0.49E- 


01 


0.065 


0.51E 


03 


PI 024- 2 


1024 


0.49E- 


01 


0.065 


0.51E 


03 


P1024-3 


1024 


0.49E- 


01 


0.065 


0.51E 


03 


P2048-1 


2048 


0.39E- 


01 


0.082 


0.26E 


03 


P2048-2 


2048 


0.39E- 


01 


0.082 


0.26E 


03 




2048 


39E- 


01 


0.082 


0.26E 


03 


P2048-4 


2048 


0.39E- 


01 


0.082 


0.26E 


03 




2048 


0.39E- 


01 


0.082 


0.26E 


03 


P4096-1 


4096 


0.25E- 


01 


0.10 


0.13E 


03 




4096 


25E- 


01 


0.10 


0.13E 


03 


P40Qfi-S 


4096 


0.25E- 


01 


0.10 


0.13E 


03 


P8192-1 


8192 


0.25E- 


01 


0.13 


0.64E 


04 


P8192-2 


8192 


0.25E- 


01 


0.13 


0.64E 


04 


P8192-3 


8192 


0.25E- 


01 


0.13 


0.64E 


04 


P16384 


16384 


0.20E- 


01 


0.16 


0.32E 


04 


P32768-1 


32768 


0.16E- 


01 


0.2 


0.16E 


04 


P32768-2 


32768 


0.16E- 


01 


0.2 


0.16E 


04 


P32768-3 


32768 


0.16E- 


01 


0.2 


0.16E 


04 


P32768-4 


32768 


0.16E- 


01 


0.2 


0.16E 


04 


P32768-5 


32768 


0.16E- 


01 


0.2 


0.16E 


04 


P65536 


65536 


0.12E- 


01 


0.25 


0.80E 


04 


P131072 


131072 


0.99E- 


02 


0.33 


0.40E 


05 


P262144 


262144 


0.78E- 


02 


0.41 


0.20E 


05 



Table 1. Details of the simulations: N is the number of randomly 
distributed particles in a sphere of diameter taken equal to unity, 
£ is the average interparticle separation (see text for definition), e 
is the softening length in the gravitational potential and m is the 
particle mass, in units in which the mean mass density is unity. 



modified at smaller scales to give a force which is attractive 
everywhere and vanishing at zero separatioijf]. 

The values of the ratio e/£ given in Table [1] correspond 
in fact to the single fixed value e = 0.0028, i.e. the smooth- 
ing in this set of simulations is fixed in units of the ini- 
tial size of the the system. Given that our aim here is to 
reproduce as closely as possible the evolution of the parti- 
cle system (i.e. without smoothing), we will test carefully 
for the dependence of our results on the choice of e, and, 
specifically, for their stability when e is decreased compared 
to the value chosen in this set of simulations. As this dis- 
cussion is closely related to the issue of the validity of the 
mean field (Vlasov Poisson) approximation to the dynamics 
of the system, we will present these results in Sect. [5] where 
we discuss this question. For now we note simply that if 
such a mean field approximation is valid, one expects that 
it should be sufficient that e be significantly smaller at all 
times than the length scales relevant to this dynamics. We 
will see that the collapse phase leads to a minimal length 
scale of the structure of orde r the scale £ (and, as estab- 
lished bv lAarseth et al.1 (119881 ) and lBoilv etai] (|2002l '). pro- 
portional to this scale). Given the values for e/£ shown in 
Table [T] it follows that the chosen e is, in the largest TV 
simulation, still smaller than this minimal collapse scale. 

The large compression which the system undergoes as it 



^ In these units therefore G = 37r/32. In physical units Tgcm is 
approximately 2100 seconds if po is Ig/cm'^. 



^ The exact expression for the smoothing function may be found 
in lSpringel et al.l 1I2OOII) . 
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Figure 2. The virial ratio for the particles with negative energy 
[see Eg JlSI l] as a function of time for different simulations. In the 

Figure 1. Total energy in the indicated simulations as a function jj^gg^ panel is shown a zoom of the behaviour after the collapse. 

of time. 



3.3 Basic qualitative and quantitative results 



collapses requires that particular care be taken in choosing 
the numerical parameters which control the precision on the 
force calculation and time stepping. Given the physics of the 
system, we have used prescriptions for parameter choices in 
which we divide the entire range of times in three different 
phases. During the first (0 t ^ 0.95Tscm) and the third 
(1.14Tscm !^ t =5! 3.5Tscm) wc havc used a relatively large 
time step (of the order of 5 x 10~''rscm) while in the col- 



lapse phase (0.95rs, 



^ t < I.Uts, 



sec) we have used 



a more accurate time step (of the order of 5 x 10~ Tscm)- 
We have performed several tests to check the stability of 
the results given below, at the relevant level of numerical 
precision. We have in all cases quantified carefully the con- 
servation of the total energy, which experience has shown to 
be a very sensi tive quantity for monitoring the accuracy of 
the simulation l|Aarsethll2003l ). Shown, for example, in Fig.[T] 
is the temporal evolution of the total energy Et normalised 
to the initial energy Eq in three of the simulations in Ta- 
ble [1] For larger the energy is conserved less well, but 
the fluctuations are of less than one percent for all but the 
largest A''. For the very largest A*" simulation in Table [T] the 
fractional error in the total energy reaches about 5% during 
the collapse phase. However, in this case and all the simula- 
tions, these variations in the energy are fluctuations about 
the total energy which do not lead to any measurable drift in 
the average value, which would be indicative of systematic 
error which propagates in timeQ. 

A further test of our simulations is their consistency 
with previous results reported in the literature for the same 
system, which we will discuss in the next subsection. We 
will then present our results for the dependence on initial 
conditions (i.e. on A'') of the final virialized state. 



The leap-frog time-integration method used by GADGET2 has, 
due to its simplectic nature, the property that it conserves ex- 
tremely well — despite local fluctuations in tim e — the t o tal en - 
ergy of the system. For a detailed discussion see ISpringell ||2005|) . 



Qualitatively the behaviour of this specific kind of cold col- 
lapse initia l conditions is well k now n, and has b e en stu died 
in detafl in lAarseth et al.l (|l988f ) and lBoilv et al.l (|2002f ): the 
system collapses following approximately the spherical col- 
lapse solution until it reaches some minimal size, at a time 
of order Tscm- It then "turns around" and rapidly (i.e. in 
a time considerably shorter than Tscm) settles down to a 
quasi-stationary virialized state. This is illustrated in Fig. [2] 
which shows, for the simulations indicated (cf. Table [T|, the 
evolution of the ratio 

^(*) = (15) 

where K" is the kinetic energy of the particles with negative 
energy, and W" the potential energy associated with the 
same particles. The ratio represents the virial ratio of the 
entire system at early times, and that of the bound particles 
at late times. 

A typical evolution of the radial density profile is shown 
in Fig.O until close to the maximal collapsed configuration, 
it maintains, approximately, the top-hat form of the original 
configuration and then in a very short time changes and 
stabilizes to its asymptotic form. We will discuss below in 
further detail this latter form and will also study in detail 
the deviations from the original form during the collapse 
phase. 

Studies of the A'^ dependence of the evolution have fo- 
cussed mainly on the question of how the minimal size 
reached by the structure depends on this parameter. This 
minimal radius Rmin may be defined in different ways, e.g., 
as the minimal value reached by the radius, measured from 
the centre of mass, enclosing 90% of the mass. Alternatively 
it can be estimated as the radius inferred from the potential 
energy of the particles, the minimal radius corresponding to 
the maximal negative potential energjjf]. The behaviour of 
Rmin, determined by the flrst method, as a function of A' is 
shown in Fig. |4l This observed behaviour is in good agree- 
ment with the prediction of the scaling Eq. (|12|) . in accord 



° Spherical symmetry is maintained up to a very good approxi- 
mation until this time. 



Energy ejection in cold spherical collapse 7 




Figure 3. Density profile of particles with negative total energy 
at different times during the collapse, for the P65536 simulation. 




Figure 5. The absolute value of the minimal value reached by 
the potential energy (i.e. at the time t = tmax), as a function of 
N. 




Figure 4. Behaviour of the minimal radius Rmin attained, de- 
termined as described in text, as a function of N. The solid line 
is the best fit to the prediction of Eq. II12I I. 




Figure 6. The behaviour of /^(t), the fraction of the particles 
with positive energy, as a function of time for two different sim- 
ulations. A dependence on the number of particles is manifest. 



with t he findings of lAarseth" et al.l (|l988l ) and iBoily et al.l 
l|2002t j^. We note that the factor of proportionahty be- 
tween Rmin and £ is in fact almost unity (as in our units 
£ = 0.5/N^^^). 

In Fig. [5] the absolute value |Wmin| of the potential 
energy at its minimum value, as a function of A'^. We again 
observe an excellent fit to the predicted behaviour. 

3.4 A'^ dependence of ejected mass 

Our focus in this paper is on the presence of an ejected com- 
ponent of the mass, and energy, which has not been closely 
examined in these previous works. One of the features of 
the collapse and virialization is indeed that while all par- 
ticles start with a negative energy, a finite fraction end up 
with a positive energy. Given that they move, from very 
shortly after the collapse, in the essentially time indepen- 
dent potential of the virialized (negative energy) particles, 
they escape from the system. Indeed at the times at which 

^ The latter paper finds that this result extend to A' of order 10^, 
an order of magnitude beyond the largest N we report here. 



we end our simulations all such particles are very far from 
the collapsed region and move outward with almost constant 
energy, with a negligible probability of an encounter which 
could stop their escape to arbitrarily large distances. 

This transfer of energy, which we will examine in detail 
below, occurs in a very short time around the maximal col- 
lapse, and depends on A'^. Both these facts can be seen in 
Fig. [S] which shows a plot of the fraction of the particles 
with positive energy at any time. 

Note that, while our simulations indicate that f at- 
tains a well defined asymptotic value, this is not expected 
to be the truly asymptotic value of this quantity, or of any of 
the other quantities we consider below: on much longer time 
scales than those considered here (i.e. the times characteris- 
tic of the violent relaxation of the system) further particles 
may gain energy and esca pe from the system, notably by two 
body encounters (see e.g. lBinnev fc Tremaind ([199J)). If the 
system is ergodic, simple considerations based o n the micro- 
canonical entropy (see e.g. lPadmanabhanI (|l990l )) imply that 
at asymptotically long times, and for an unsmoothed gravi- 
tational potential, the particles will tend to a configuration 
in which there is a single pair of particles with arbitrar- 



8 M. Joyce, B. Marcos and F. Sylos Labini 



y = 0.048 + 0.022 ln{x) 



A' 

Figure 7. Behaviour of the fraction of ejected particles as a func- 
tion of the total number of particles in the system. The solid line 
is the phenomenological fit given by Eq. I|16| l. 



ily small separation, and the rest o f the ma s s is in an ever 
hotter gas of free particles (see e.g. lAarsethl ()l974h ). When 
the potential is smoothened, as in our simulations here, we 
would expect on such very long time scales to obtain a final 
state with some part of the mass bound and which depends 
strongly on this smoothing. We will, however, not explore 
this temporal regime here. 

The A'^ dependence of the asymptotic value of we 
identify in this way is shown in Fig. [7] We note that, al- 
though determinations in different realizations of a given 
system with the same number of particles fluctuate, there 
appears to be a very slow, but systematic, increase as a 
function of A*', which is reasonably well fit by 



/f(iV) f=ia + blog{N) 



(16) 



where a = 0.048 and b = 0.022. Alternatively it can be fit 
quite well (in the same range) by a power law f « O.IN^'^ . 

3.5 A'^ dependence of ejected energy 

While the factor of the mass ejected varies, at most, very 
slowly with A'', the energy it carries away has a much 
stronger dependence on A'^. In Fig. |S] is shown, for exam- 
ple, the behaviour of the total kinetic energy as a function 
of time in two simulations with different A'^. Not only the 
time at which the maximum is reached, but also the final 
value of the energy, change clearl\F°l. We consider now care- 
fully this dependence, and the related one of the potential 
energy of the different components. 

In general the total energy, which is equal to the initial 
energy Eq, may be written as 

Eo^W + K = W'' + W'' + K'' + K'' + W^" (17) 

where W {W, W") and K {K", K") are the total po- 
tential and kinetic energy of the particles (with positive 
energy, with negative energy) at the time t, and W^^" is 
the potential energy associated with the interaction of the 
particles with positive energy with those of negative total 




Figure 8. Behaviour of the total kinetic energy for two simula- 
tions with different number of particles. 




Figure 9. Upper panel: behaviour of the kinetic and potential 
energy, normalized in units of the absolute value of the initial 
energy, for all particles and for particles with negative total energy 
as a function of time for the P65536 simulation. Lower panel: The 
same for particles with positive total energy. 



energy. Note that Eq is simply, up to a correction from 
fluctuations due to the particle discreteness, the gravita- 
tional potential energy of a uniform ball of radius ro, i.e., 
Eo = W{t = Q) = -f where M = mN is the mass of 

the ball (see e.g. iBinnevfc Tremaind (|l994l )). 

Shown in Fig. |9] are the evolution of these guantitief^. 
The maximum of the kinetic energy K corresponds evidently 
to the minimum of the potential energy W. We see also that 
W ~ W" , with only a very slight deviation around the time 
of collapse. This corresponds simply to the expulsion of the 
positive energy particles, which makes both and W^^" 
completely negligible once these particles are far from the 
remaining bound part of the system. 

Once the collapse phase is over, we therefore have to an 
excellent approximation 



(18) 



^'^ The common evolution of both curves is simply that predicted The energies here and in subsequent figures, unless specified, 

by the spherical collapse model, which leads to a divergence of the are given in units of the absolute value of the initial total energy 

kinetic energy at Tscm- Eo- 
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Figure 10. Behaviour of the total kinetic energy of the ejected 
particles (normalized to initial total energy) as a function of the 
number of particles in the system. The dashed line is the best 
fit to a single power-law behaviour for RP itself, the solid line a 
fit to Eq. (|20) with W" oc N'^^^, i.e. with W" scaling exactly as 
^min observed above. 



Further, since the bound particles rapidly virialize into a 
quasi-stationary state, we have 



2K" + W" ^0 . 



(19) 



The three non-zero energies K^, and can there- 
fore be expressed in terms of the initial (potential) energy 
Eq and a single a priori unknown quantity, e.g., the ejected 
energy . 

Note that it follows from Eq. gill and Eq. ^ that 



W" 



(20) 



so that for the case we are studying, in which Eq < 0, these 
considerations do not on their own imply the necessity for 
energy ejection, i.e., one can satisfy Eq. (|2U|I with = 0. 
For Eo > 0, however, the formation of a bound virialized 
stationary state requires such ejection. 

In Fig. [To] we show the observed behaviour in our sim- 
ulations of the total kinetic energy of the ejected par- 
ticles as a function of N. Two fits are shown: one is to a 
single power-law behaviour for K'' itself, the other to the 
behaviour which would be observed, using Eq. (|20p . if W" 
scaled as VP""'" observed above. Thus in these curves we see 
a behaviour which is close to, but clearly different from, the 
simple scaling observed for VF"**". 

Shown in Fig. [Tl]is instead the ratio K^/f^, i.e., the 
kinetic energy per unit ejected mass, as a function of N. We 
observe that the simple behaviour / oc A'^^^^ provides 
a very good fit, with less scatter than in the previous plot. 
After we have discussed below in detail the mechanism of 
mass/energy ejection, we will give a simple scaling argument 
which accounts for this behaviour. 



3.6 A'^ dependence of density profiles 

The radial density profile of the virialized structure formed 
by the bound particles after the collapse may be well fit in 
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Figure 11. Observed behaviour of the ratio / fP for the set of 
simulations in Table [T] 
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Figure 12. Behaviour of the parameter ro as a function of the 
number of particles in the simulation. 



all our simulations by the functional form 

, s no 
n{r) 

'l+ 



(21) 



While the previous studies of exactly cold uniform initial 
conditions do not report results for this quantity, we ob- 
serve that it is in agreement with that found for cold (i.e. low 
initial virial ratio) initial c onditions in several previous stud- 
ies (se e e.g. iHenonI l| 19641 ): Ivan Albadal (|l982l ): IRov fc Perez! 
(|2004l )'). 

While the form of this profile is observed to be very 
stable in our set of simulations, the parameters no and ro 
vary as shown in Figs. ll2l and [T3l Thus a good fit is given by 
the simple behaviours ro oc A'^^^''^ and no oc A'^^. In Fig. [TJ] 
we show the density profiles for various simulations with 
different A'^ where the axes have been rescaled using these 
behaviours. 

These behaviours can be easily related to those we have 
observed in the previous section for and . Indeed us- 
ing the ansatz of Eq. (1211) it is straighforward to calculate 
both the number of particles TV" = (1 ~ f^)N which are 
bound, and the binding energy W" of these particles [which 
is then related to through Eq. ([20}]. The first may be 
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Figure 13. Behaviour of the parameter no as a function of the 
number of particles in the simulation. 
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Figure 14. Density profile of the virialized structure at a time 
t 4:Tscm for simulations with different number of particles. The 
y-axis has been normalized by A'^^ and the x-axis by N~^^^ (see 
text for explanations). The behaviour of Eq. II21I I is shown for 
comparison. 



determined analytically as 
while the second may be written as 



AGm^nnrn 



(22) 



(23) 



where ~ 11 is a numerically determined constanlF^. 

The fitted behaviours for ro and no thus correspond, 
since m oc 1/N, to A'^" ~ A'^ (i.e. a constant bound mass, 
and therefore a constant ejected fraction of the mass f^) 
and W" ~ N^^^ (and therefore oc N^'^). Modifications 
of these fits consistent with the very slow variation of de- 
scribed in the previous subsection can be given. Thus, to a 
good approximation, we simply find that the characteristic 
size of the final structure scales with A'^ as Rmin, the mini- 
mal radius attained in the collapse, does, i.e., in accordance 



N 


ti 


t2 




t4 


131072 


0.904 


0.988 


0.998 


0.999 


8192 


0.666 


0.952 


0.988 


0.998 


512 




0.809 


0.953 


0.993 



Table 2. Times chosen for the different number of particles. 



with the simple scaling argument which gives the predic- 
tion Eq. p2p . Correspondingly the potential energy which 
is bound in this structure, given that the mass is approx- 
imately constant, increases in proportion to the inverse of 
this scale (and scales also as minimal value of the potential 
energy reached during the collapse). 



4 MECHANISM OF MASS AND ENERGY 
EJECTION 

In this section we consider in detail how the mass is ejected 
from the collapsing system. 



4.1 The ejected particles 

An evident starting point in considering the ejection of mass 
is to ask whether there is a direct correlation between a par- 
ticle's initial radial position and the likelihood of its ejection, 
i.e., are particles preferentially ejected from the outside or 
the inside of the initial sphere? 

Shown in Fig. [15] is the "radial rank correlation" 0: at 
the indicated time t one ascribes a rank n{t) € [1, A'^] to each 
of the A*' particles according to their radial distance from 
the centre of mass, and then for each point one plots the 
couple {n{t)/N, n{0)/N). The chosen times ti...t4 are not 
the same for the different N, but are as given in Tabled This 
choice of the times is that discussed in Sect. l2.3l above. i.e., it 
corresponds to comparing the different simulations at times 
at which they are equivalent in the approximation — which 
leads, as explained in Sect. 12.31 to the predicted observed 
scaling of Rmin — that they behave like infinite continuous 
systems with fiuctations initially of amplitude proportional 
to l/\/iV, i.e., the times t' for the particle number N' > N 
are determined by 



(24) 



where f{R) is the growth factor for perturbations, of which 
the full expression is given in Eq. (|A3|I . 

From these plots we see that that, while the strong cor- 
relation evidently present at early times tends to disappear, 
there is still some visible structure in the plot even at the 
final time, which corresponds to a scale factor just slightly 
larger than Rmin, i.e., just a very short time before the max- 
imal collapse. We do not show the plots for later times — 
after the collapse — as they are essentially identical to this 
last one. Given that typically of order twenty per cent of the 



^■^ For an isolated spherical system with a mass density profile 
p(r) the potential energy is —AnG Jq°° M{r)p{r)rdr, where M{r) 
is the mass enclosed in the radius r. 



Th is statistic has been introduced and used bv lAarseth et all 
lll988l) to examine the degree of correlation of particles' radial 
positions with their initial radial positions. From their data for 
N = 5000 it is concluded that such correlation is very weak. 
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Figure 16. Normalized liistogram (i.e. estimated probability dis- 
tribution function) of the ejected particles as a function of their 
radial position in the initial configuration. The dashed line shows 
the best power-law fit. 

mass is ejected, it thus appears that although ejected mass 
originates from all parts of the initial structure, it comes very 
preferentially from the outer regions of the sphere. Indeed 
the residual correlation in the plot appears to be a result 
solely of the correlation between ejection and initial radial 
position. 

The degree of this correlation can be seen much more 
clearly in Fig. 1161 which shows the estimated probability 
density function of ejection as a function of initial radius, 
i.e., the normalized histogram of the ejected particles as a 
function of their initial radius. In all the simulations we ob- 
serve a very consistent behaviour, fitted approximately by a 
simple power-law: 

Pe(r) oc . (25) 

Given the uniform initial distribution, the probability that 
a particle is at radius r initially is proportional to r^, so this 
result corresponds to a conditional probability of ejection 

p{r) oc , (26) 

i.e., the probability that a randomly chosen particle at ra- 
dius r in the initial configuration will be ejected grows ap- 
proximately as r^. Thus, although some of the ejected mass 
comes from the inner parts of the initial sphere, there a very 
clear systematic correlation between initial radial position 
and ejection. 

4.2 Development of lag during collapse 

Since particles are preferentially ejected from the initially 
outer parts of the sphere, one would expect that their ejec- 
tion may be related to a dependence of the evolution of par- 
ticles prior to the collapse on their radial position. Indeed we 
have noted that the radial rank correlation plot of Fig. [15] 
at the last time shown, just before the collapse, is almost 
indistinguishable from that after the collapse. As we have 
discussed, however, in the uniform limit of the spherical col- 
lapse model (SCM), "particles" do not "see" the finite size 
of the system during the collapse phase, and indeed pertur- 
bations about the SCM in the continuum (fluid) limit may 
be treated in this same approximation that the collapsing 



system is infinite. With these approximations therefore the 
evolution of the trajectory of a particle cannot depend non- 
trivially on its initial radial position. As a corroUary such 
a dependence on radial position must arise necessarily from 
a coupling between the evolution of perturbations and the 
finite size of the system. 

To quantify the difference in evolution of particles as a 
function of their initial radial position we consider the "lag" , 
which we estimate as follows: 

^ ('■'W-^scmW), (27) 

7'^ I r' I ^ r + (5r 

where the sum is over the particles N{r) initially at radial 
distance between r and r + Sr (with respect to the centre of 
mass), r'{t) is actual radial distance at time t of the particle 
and r'gQyi{t) is the radial distance predicted at time t in the 
SCM, i.e., rscM(i) = R{t)r'{0). The lag thus represents the 
average discrepancy between a particle's radial position and 
that in the SCM, as a function of radius. 

The results for this quantity for difi'erent simulations 
and times — the same as those in the previous figure of 
the radial rank correlation — are shown in Fig. 1171 For 
convenience of comparison A{r, t) has been normalized to 
the radius rscM(i) of the sphere in the SCM model at the 
corresponding time t. We see that a net positive lag (i.e. 
net "stretching" of the positions with respect to their SCM 
values) develops initially for two quite separate ranges: for 
particles initially close to the centre, and for particles close 
to the outer boundarjj^. The reason for the development of 
systematically positive values is quite different in the two 
cases. For the particles close to the centre there is a sys- 
tematic effect which makes the estimator necessarily biased 
towards positive values: once particles can deviate from their 
SCM radius by of order Ar, particles within this distance 
of the centre can reach the centre of mass, i.e., the minimal 
possible value of the radius, and thus their radius will start 
to increase again, which leads to an intrinsic asymmetry to- 
wards net positive values of A(r, t) at these scales. From the 
corresponding plot in Fig. [15] we see that particles' radial po- 
sitions with respect to the SCM indeed grow in the course of 
the collapse, with considerable redistribution of mass over 
the entire system already at the third time slice. 

The clear positive lag which develops first close to the 
outer boundary, and then propagates progressively into the 
volume as time goes on, can be understood as follows. In 
addition to the mean field responsable for the SCM mo- 
tion, the particles move under the effect of the field due to 
fiuctuations which modify the SCM trajectories. These lat- 
ter contributions are local, with the dominant contribution 
coming initially from nearest neighbour particles. For a shell 
at the outside of the volume, there is thus a difference com- 
pared to one in the bulk: as particles move around there is 
no compensating inward fiux at the boundary for the mass 
which moves out under the effect of perturbations. Thus the 
net density of the outer shell decreases, and also the average 
density in the sphere at the corresponding radius, slowing 



Note that, since the mass inside the radius r grows in propor- 
tion to r'^, the fraction of the mass lagging in the two ranges is 
actually comparable even at the initial time despite the difference 
in radial range. 
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Figure 15. Radial initial rank (horizontal axis) versus rank at the times indicated in Table [2] for, from left to right, N = 131072, 
N = 8192 and TV = 512 particles and, from top to bottom, the times ti, t2, <3 and t^. For the plots corresponding of the system with 
N = 131072 particles we have performed a random selection of 8076 particles. 



its fall towards the origiiQ- As time goes on this asymme- 
try propagates into the volume: because the outer shells fall 
inward more slowly the mass in these shells stretches out 
radially, lowering the density further, which in turn "feeds" 
less flux into the shell. In the subsequent times in the figure 
we see clearly this effect, with a larger part of the outer mass 
lagging more and more compared to the SCM model, and, 
more significantly, compared to the mass closer to the cen- 
tre. Indeed at the last time shown, just before the minimal 
collapse radius is reached, the lag is essentially constant — 

We note again, as we did in the introduction, that the forma- 
tion of a lagging low density outer region ( "halo" ) during collaspc 
is noted in Aarsoth et al. (1988), and a similar qualitative expla- 
nation for it is briefly outlined. 



equivalent to a uniform dilatation of the whole sphere with 
respect to the SCM — for particles initially at radii less than 
approximately 0.8, but then rises sharply. The conclusion is 
that particles in the outer shells arrive at the centre of mass 
on average much later than those in the bulk. 

4.3 Mechanism of mass ejection 

That it is this "late arrival" of the outer parts of the sphere 
which leads to their net gain in energy and subsequent ejec- 
tion can be seen from Fig. 1181 It shows, for the P131072 
simulation, the temporal evolution of the components of the 
mass which are asymptotically ejected or bound (i.e. with 
positive or negative energy a short time after the collapse 
has occurred). More specifically it shows the evolution of «e 
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Figure 17. The "lag" as defined in Eq. I|27| l. normalized as indicated, for the same four times as in tlie previous figure. 



(and Vb) which is the average of the radial component of the 
velocity for the ejected (and bound) particles, and also Ce 
(and et) which is the mean energy per ejected (and bound) 
particle (i.e. the average of the individual particle energies). 
The behaviours of Ve and Vb show clearly that the ejected 
particles are those which arrive on average late at the cen- 
tre of mass, with Ve reaching its minimum after the bound 
particles have started moving outward. Considering the en- 
ergies we see that it is in this short time, in which the for- 
mer particles pass through the latter, that they pick up the 
additional energy which leads to their ejection. Indeed the 
increase of Be sets in just after the change in sign of Vb, i.e., 
when the bound component has (on average) just "turned 
around" and started moving outward again. The mechanism 
of the gain of energy leading to ejection is simply that the 
outer particles, arriving later on average, move through the 
time dependent decreasing mean field potential produced by 
the re-expanding inner mass. We note that a mechanism of 
th is kind for "mean - field" mas s ejection has been discu ssed 
in lDavid fc Theund (|l989l ) and lTheuns fc DavidI (|l990l ). Al- 
though we cannot quantitatively apply the results of these 
latter works to the present case — they treat the case of 
oscillations about a quasi-equilibrium — they give a quali- 
tative insight into the ejection we observe. 
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t 
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Figure 18. Radial velocity, and average energy per particle, as a 
function of time, of particles which are bound/ejected at the end 
of a P131072 simulation. The energy of the particles has been 
arbitrarily rescaled. 



4.4 Scaling of ejected energy 

Let us now finally attempt to be a little more quantitative, 
and relate this qualitative understanding of how the mass is 
ejected to the actual scaling we have observed of the ejected 
energy as a function of A'^. 

To a first approximation the curves at the latest time 
plotted in Fig. [17] overlap, i.e., the profile describing the 
lag of the outer shells with respect to the inner ones is the 
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Figure 19. Time interval of the "peak" of the kinetic energy as 
a function of A'^. 




Figure 20. Time t = tmax — Tscm as a function of A'. 



same. In this approximation all the simulations are therefore 
equivalent, up to rescalings of time and length, and would 
be expected to eject the same fraction of their mass if this 
ejection arises from the delay described by these curves. 

A simple estimate of the energy Eg gained by the outer 
lagging mass as it passes through the inner core is then 



Rn 



tc 



'^char 



(28) 



where the first term represents the rate of change of potential 
produced by the "core" structure, and the second the typical 
time for particles to cross it, where the parameters t^har 
and Vchar are characteristic time and velocity scale for the 
process. We have seen already numerically that, to a very 
good approximation, W, 
simply 



oc R^i„, and therefore we have 



Ea 



tchar 



(29) 



We expect the role of tchar to be played by the typical 
time over which the "fall through" leading to the transfer of 
energy between the components occurs as we have described 
above. This can be estimated numerically, and its scaling 
with A'^ inferred, from the data we have shown. For example, 
estimating it from the curve for the evolution of the total 
kinetic energy, which reaches a sharp peak as we go through 
this phase, as the temporal width at half the maximal value 
of the kinetic energy, we find the results shown in Fig. 1191 
As shown a very good fit is given by the scaling At oc N^^^^ 

The scaling of this characteristic time coincides with 
that predicted by the perturbed SCM model, as described 
in Sect. 12.31 for the difference tmax — fscm, where tmax is the 
time at which the maximal collapse is reached. That this 
latter quantity also scales itself as predicted can be seen in 
Fig. [20] 

The characteristic velocity v^har is that of the particles 
which will be ejected as they pass through the collapse, i.e., 
Ve in Fig. [18] above. We might anticipate that its scaling 
will also be predicted by the SCM for the maximal velocity 
reached at the collapse, [cf. Eq. (|14|l above], i.e., 

Vchar OC N^"" . (30) 

That this is indeed a good approximation to the scaling can 
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Figure 21. Radial velocity, as a function of time, of particles 
which are bound (solid lines) or expelled (dashed lines) at the end 
of the different simulations indicated. For the 8192 and 131072 
simulations, the x-axis and y-axis have been rescaled as described 
in the text. 



be seen in Fig. 1211 which shows the same two radial velocities 
He and «(, as in the previous figure. Both axes of the plots for 
the curves from the P512 and P8192 simulations have been 
rescaled: the time axis as just described, and the velocities 
according to Eq. (|30p . We see that in these rescaled variables 
the evolution is indeed very similar. 

Taking these scalings to be exact, Eq. (|29|) then gives 



Ea^N 



1/3 



(31) 



This simple estimate agrees quite well with the observed 
scaling of the ejected energy (i.e. asymptotic kinetic energy 
of the positive energy particles) . We have seen, however, 
that this scaling is followed more precisely by / , i.e., 
the energy ejected per unit mass, while has a very weak 
(approximately logarithmic) dependence on A''. The simple 
result obtained here has in fact been obtained with the as- 
sumption that the mass ejected is independent of A'^. As- 
suming the relevant properties for understanding the ejected 
energy are given by the curves for the lag shown in Fig. 1171 
this corresponded to the approximation that the curves at 
the latest time, just before the collapse, are the same for all 
A'^. In fact we can see that this is only true to a first approx- 
imation, and specifically that the integrated "lagging" mass 
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(i.e. under the curve ) indeed appears to grow slowly as a 
function of A'^. Thus the scaling result obtained above can 
be adapted taking into account this increase as a function 
of A'^ of the mass. We do not, however, have an analytical 
understanding of this observed A'^ dependence of the lagging 
(and subsequently ejected) mass. 

5 THE VLASOV POISSON LIMIT AND 

SENSITIVITY TO INITIAL CONDITIONS 

5.1 Definition of the Vlasov Poisson limit 

We now discuss the extent to which the evolution of the sim- 
ulated system is "collisionless" , i.e., described by the coupled 
Vlasov-Poisson (VP) equations (or "collisionless Boltzmann 
equation"). The latter should describe the evolution of the 
system in an appropriate N oo. The naive such limit 
applied here, i.e. N ^ oo Poisson distributed particles, is 
clearly not the desired limit, as it converges to the singu- 
lar SCM. Indeed we have explicitly identified macroscopic 
A*' dependences in various quantities, which diverge if we 
perform this naive extrapolation. 

Existing formal p roofs of the validity of the VP limit 
l|Braun fc Hepdll977l ) for a self-gravitating system require, 
however, that the singularity in the gravitational force at 
zero separation be regulated when the limit A?^ — > oo is 
taken Incorporating such a modification of the force will 
clearly regulate the divergence in the SCM as we converge 
to the uniform mass distribution. One would then expect to 
obtain, for sufficiently large A', a final state which is well 
defined and A'^ independent, but strongly dependent on the 
implementation (and scale) of the regulation. 

This limit also is not the VP limit relevant here: while 
we have indeed introduced such a regulation of the force 
(characterized by the smoothing e), we have done so, as 
discussed in Sect. |3] for reasons of numerical convenience. 
Indeed, as we have discussed, our criterion for our choice 
of £ is that it be sufficiently small so that our numerical 
results are independent of it, and we interpret our results 
as being representative of the limit e = 0. In Fig. [21] are 
shown, for example, the evolution of the fraction of particles 
with positive energy as a function of time for the different 
indicated values of e. Other quantities we have considered 
show equally good convergence as e decreases. Note that for 
the given simulation (cf. Table [T]), the mean interparticle 
distance £ = 0.016, so that the convergence of results is 
attained once e is significantly less than £. As we have seen, 
the minimal size reached by the collapsing system scales 
in proportion to £ (with a numerical factor very close to 
unity, as can be seen in Fig. [4]). We interpret the observed 
convergence as due to the fact that the evolution of the 
system is determined primarily by fluctuations on length 
scales between this scale and the size of the system. Once e 



° This proof is given in a finite volume. We do not believe that 
this is an important difference here, as the dynamics we describe 
will be modified in a quite trivial way, on the time scales con- 
sidered, if we put our system in a box: the ejected component 
will simply bounce off the walls and remain as an unbound cloud 
moving in and out of the time independent potential of the "core" . 




Figure 22. Evolution of the fraction of the mass with positive 
energy for simulations with N = 32768 for the different values 
indicated of the smoothing parameter e. 

is sufficiently small to resolve these length scales at all times, 
convergence is obtained. 

A different, less rigorous but more physically in- 
tuitive, derivation of the VP limit is given (see e.g. 
iBuchert fc Dominguej (j2005)) by a coarse-graining of the 
exact one particle distribution function over a window in 
phase space. The VP equations are obtained for the coarse- 
grained phase space density when terms describing pertur- 
bations in velocity and force below the scale of the coarse- 
graining are neglected. A system is thus well described by 
this continuum VP limit if the effects of fiuctuations below 
some sufficiently small scale play no role in the evolution. 
The validity of such an assumption for our system is indi- 
cated by precisely the kind of behaviour we have just noted 
of our results as a function of e. 

5.2 Numerical extrapolation to the VP limit 

The VP limit defined in this way thus corresponds to taking 
N ^ oo, while keeping fixed the initial fluctuations above 
some scale. The validity of this limit for our system may be 
explored numerically by defining such an extrapolation and 
studying stability of our results to it. This can be done as 
follows: starting from a given Poissonian initial condition of 
TV particles in a sphere of radius R, we create a configura- 
tion with A'^' = nN particles by substituting each particle 
by n particles in a cube of side 2rs centred on the origi- 
nal particle. The latter particles are distributed randomly 
in the cube, with the additional constraint that their centre 
of mass is located at the centre of the cube (i.e. the centre of 
mass is conserved by the "splitting"). In this new point dis- 
tribution fluctuations on scales larger than are essentially 
unchanged compared to those in the original distribution, 
while fluctuations at and below this scale are modifiecO- We 
have performed this experiment for a Poisson initial condi- 
tion with N = 4096 particles, splitting each particle into 
eight (n = 8) to obtain an initial condition with N' — 32768 

See lGabrielli &: Joyce! l(20o3) for a detailed study of how fluc- 
tuations are modified by such "cloud processes" . 
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Figure 23. Evolution of the fraction of particles with positive 
energy as a function of time, for the different indicated values 
of the parameter described in text. The "original" initial con- 
ditions has 4096 particles while the others have 32768 particles. 
Note: the curve for 0.8£ is not visible because it is superimposed 
on that for the "original" one . 
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Figure 24. Evolution of the fraction of particles with positive 
energy as a function of time, from initial conditions given by 
the "split" configurations described in text. The results are for a 
Poisson configuration with = 4096 split into one with N' = 
32768 particles, and = 0.8£, and differ only in whether the 
centre of mass of the split particles is conserved or not. 



particles. Results are shown in Fig. [23] for the ejected mass 
as a function of time, for a range of values of the parameter 
Vs , expressed in terms of £, the mean interparticle separation 
(in the original distribution). While for Vs = 0.8£ the curve 
of ejected particles is indistinguishable in the plot from the 
one for the original distribution, differences can be seen for 
the other values, greater discrepancy becoming evident as 
Ta increases. This behaviour is clearly consistent with the 
conjecture that the macroscopic evolution of the system de- 
pends only on initial fluctuations above some scale, and that 
this scale is of order the initial interparticle separation £. 
And, as anticipated, this translates into an A'' independence 
of the results when N is extrapolated in this way for an Vs 
smaller than this scale. 

The observed behaviour of the ejected mass for the 
larger values of can be understood as follows: for an ini- 
tial Poisson distribution, after dividing each particles follow- 
ing the procedure described above, the distribution remains 
Poissonian at large scale. Therefore the growth of fluctu- 
ations inside the sphere and hence the radius of minimal 
collapse, etc., will be the same for the original and the new 
simulation. The relevant difference between the initial con- 
ditions of the two simulations are for radius r > R — Vs in the 
original distribution, because particles situated in this region 
have a strong probability to lie outside the original sphere 
in the split distribution. These particles will lag strongly be- 
cause they feel a weaker effective gravitational field than the 
ones inside the sphere and a large number of them will be 
ejected, explaining the strong increase of ejected particles as 
Ts increases. 



5.3 Sensitivity to initial fluctuations at large 
scales 

It is interesting to probe further the sensitivity of the results 
to the initial fluctuations over the range of scales which ap- 
pear to be relevant to the determination of the final state. In 
Fig.[24]are shown the same quantity as in the previous figure, 
for rs — 0.8£, with the only difference that the n points are 



now distributed randomly in a cube of side 2rs both with and 
without the constraint that the centre of mass be conserved. 
Without the constraint there is an additional random dis- 
placement of the mass, which induces slightly greater fluctu- 
ations at large scales, and specifically a non-zero fluctuation 
in the dipole moment of the whole mass distributiorS. We 
see that, although this corresponds to an extremely small 
modification of the fluctuations at large scales, it leads to a 
perceptible change in the final macroscopic properties. 

Finally we can probe the dependence on initial fluctua- 
tions at different scales by comparing the evolution from the 
Poissonian initial configuration, to that from initial "shuf- 
fled lattice" conflgurations with the same number of par- 
ticles. The latter are generated by placing particles on a 
perfect lattice, and then subjecting them, independently, to 
a random displacement in a cube of characteristic size Sa ■ A 
sphere containing the required number of particles is then 
extracted, with centre on a lattice point. If the scale 5a is 
or order the interparticle distance or larger, the latter dis- 
tribution has then, up to this scale, fluctuations which are 
Poissonian (and identical in amplitude to those in the initial 
Poissonian distribution), while at larger scales it has fluctu- 
ations of a much lower amplitude than in the latter distribu- 
tion. In Figs. [25] we show the evolution of the fraction of the 
ejected mass as a function of time, for the indicated values 
of the parameter 5a, for configurations with 4143 particles. 
In comparison with the Poisson configurations we have re- 
ported above with approximately the same number of par- 
ticles (4096), which eject 20 — 25 percent of their mass, the 
mass ejected is significantly larger in all cases, increasing 
monotonicaUy as 5a decreases. This is in line with what ev- 
erything we have seen above: for any 5s ^ 1 the fluctuations 

Fo r an infinite point distrib ution (see, e.g. iGabrielli et al.l 
||2004|) : iGabrielli fc Jovc'3 1I2OO8II ') the constrained case produces 
long wavelength fluctuations with a power spectrum decaying in 
amplitude in proportion to fe* (where k is wavenumber), while 
the unconstrained "shuffling" produces a power spectrum pro- 
portional to . 



Energy ejection in cold spherical collapse 17 



0.5 



1.5 



t 



5s=0 - 
5,=0.05 -- 
1=0.1 - 

5,=0.5 

Ss'1-0,- 

2.5 3 



3.5 



Figure 25. Evolution of the fraction of particles with positive 
energy as a function of time, for simulations from "shuffled lat- 
tice" initial conditions with different values of the parameter 5g. 
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Figure 26. Evolution of the fraction of particles with positive 
energy for the different indicated initial conditions (see text). In 
each case the thick line corresponds to the heavier particles (mass 
mi) and the thin line to the lighter particles (mass m2). 



at relevant scales are smaller than those in the Poisson dis- 
tribution, and their amplitude decreases as Ss does. For the 
limit case of a perfect lattice (i.e. Ss — 0), the collapse is 
the most violent. Indeed in this case the only density fluc- 
tuations regularizing the collapse are surface fluctuations 
(associated with the finite size of the system), and close to 
half (43%) of the initial mass is ejected. 

5.4 Further test for non-VP effects 

A useful test for non-VP effects is to use particles of dif- 
ferent masses in the simulations: given that the VP limit is 
essentially a mean field approximation, the trajectories of 
particles should be independent of their mass if it is valid. 
There should therefore be no segregation (up to statistical 
fluctuations) between the different species in the final config- 
uration if this limit applies. We have thus run three simula- 
tions with dispersion in particles' masses: we denote by DM1 
and DM2 two realizations of an A'' = 4096 Poissonian ini- 
tial condition, in which we have assigned, randomly, a mass 
mi = 4/3m to half of the particles, and m2 = 2 /3m to the 
other half. Finally, we denote by DM3 the same realization 
as DM2 (i.e. with particles in the same positions), but now 
with half of the particles with masses mi = 20/1 Im and 
the other half = 2/llm. In Fig. [^Hlwe show the fraction 
of ejected particles for the three simulations. Firstly we ob- 
serve that the difference between the number of light (m2) 
and heavy (mi) particles is, albeit larger in DM3 than in 
DM1 and DM2, in all three cases of order the dispersion in 
the average number ejected in different realizations (see also 
Fig- El - The difference in the number ejected of each species 
in a given simulation (with given density fluctuations) is, 
however, too large to be a purely statistical fluctuation, and 
indeed we see that there are systematically more light parti- 
cles ejected than heavy ones, with a clearly more pronounced 
effect as the mass difference increases. This small excess of 
lighter ejected particles can naturally be attributed to the 
ejection of some (small) part of the mass by two-body colli- 
sions, rather than by the mean field mechanism described in 
the previous section: such collisions modify the energy per 
unit mass, which is the relevant quantity in the mean field 



limit, in a different way according to their mass. In a colli- 
sion between a light particle and a heavy particle, notably, 
the trajectory of the former is deviated much more than that 
of the latter, and therefore it is more likely to be ejectecjf]. 



6 CONCLUSIONS AND DISCUSSION 

We have studied in this article the collapse and violent re- 
laxation of A'^ (classical Newtonian) self-gravitating parti- 
cles, initially at rest and distributed randomly in a sphere. 
We have focussed on the characterisation of the macroscopic 
properties of the quasi-stationary virialized state which re- 
sults, and in particular on the mass and energy of the bound 
and ejected components. We have identified unexpected non- 
trivial dependences of the latter on A'^, where the latter varies 
over approximately three orders of magnitude (from a few 
times lO'^ to a few times 10^). We have studied in detail 
the evolution through the collapsing phase, identifying the 
physical mechanism leading to this mass and energy ejec- 
tion. Further we have used this to give a simple analysis 
which allows one to explain the observed scaling with A*' of 
the energy per unit ejected mass (oc Afi/3). Finally we have 
clarified how one can test the hypothesis that the evolution 
is representative of the Vlasov- Poisson limit by an appro- 
priate extrapolation of the particle number. While our main 
convergence test show good stability of results, another one 
allows us to detect residual non-zero effects corresponding 
to deviation from this limit. 

One point which we emphasize is that we have not ex- 
plained, even qualitatively, the very slow (approximately 
logarithmic) growth with A'^ of the fraction of the mass 
ejected /''. We note that this result is numerically less robust 
than the result for the energy: the effect is both very small 
over the range of A^ simulated, and also shows greater fiuc- 
tuations from realization to realization than for the energy 

Such mass segregation effects have been observed also in 
lAarsethl l ll974) in the ejection by two body collisions which oc- 
curs as such systems relax on much longer time scales than those 
considered here. 
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per unit mass. Indeed in Fig. [7] we see that if we omit the two 
points with largest A'', the resuhs would be quite consistent 
with a convergence of the ejected mass to a constant value. 
As these largest A'' simulations are those we have not been 
able to test for discreteness effects using the tests we have 
discussed (which require extrapolation to still larger A''), we 
cannot exclude that the apparent continued growth might 
be due to such effects. Indeed we note the growth of the 
ejected mass with A'' cannot be consistently extrapolated by 
a logarithmic (or power law) dependence to arbitrarily large 
A'^, as it is a quantity which is bounded above (by the to- 
tal mass). The numerical results we have given, despite the 
fact that we have considered a very large number of parti- 
cles, thus do not actually resolve the asymptotically large 
A'^ dependence (or asymptotic constant value) of the ejected 
mass, or indeed that of the ejected energy (or the associated 
virialized state) is. Further numerical study, with larger sim- 
ulations, would be feasible for groups with greater numerical 
resources than ours, and might resolve the issuj^. It is nat- 
ural to expect that the slow growth in the mass will reach an 
upper bound when the fraction of mass becomes of order the 
total mass. We remark that a log A^ behaviour, for example, 
could be explained if the quantity of mass which "lags" as 
we have described, and is then ejected, grows exponentially 
in the course of the evolution from a fluctuation which is 
originally of order 1/N (or some power thereof). One would 
expect the growth law to change, and flatten to a constant, 
when the fraction of lagging mass becomes of order one half. 

In relation to this last point it is interesting also to 
consider a little further the case of an initial lattice-like con- 
figuration, which we discussed briefly in Sect. 15.31 These 
are, as we have noted, much more uniform configurations 
in the range of scales which are relevant to the macroscopic 
evolution, and indeed we observed considerably larger mass 
ejection. This corresponds, as one would expect, also to a 
more violent collapse, with larger energy ejection. We show, 
for example, in Fig. [57] the evolution of the total kinetic 
energy for the different values of 5s (the "shuffling" of the 
lattice) indicated. The increasing violence of the collapse is 
evident as 5s decreases. We note that, for sufficiently small 
5s, the maximal collapse (corresponding to the maximum of 
the kinetic energy) actually occurs before the time predicted 
by the uniform SCM. This is in contrast to the behaviour 
for all the Poissonian initial conditions we have considered, 
which in all cases gave a collapse time longer than Tscm (cf. 
Fig. I20|) . In the latter case this behaviour can be under- 
stood in terms of the lagging mass we have discussed, which 
leads to the system having a lower effective mean density 
(and therefore a longer collapse time). That a different be- 
haviour occur for the exact (or near exact) lattice, may be 
indicative that a quite different behaviour of the Poissonian 
system may be reached when A'^ is extrapolated sufficiently 
far that, like in the exact lattice, a regime is reached in which 
the dominant fluctuations are those arising from the finite 
size of the system. A naive extrapolation of the A' depen- 
dence of for the Poissonian configurations to the value 
observed for the lattice gives A^ ~ 10^. To resolve this point 
it would evidently be interesting to study the evolution from 
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Figure 27. Evolution of the kinetic energy as a function of time, 
for the shuffled lattice initial conditions with the indicated values 
of (5,, . 



these initial lattice or lattice-like configurations, and indeed 
other very un i form c onfigurations such as those described in 
[Hansen et al.l (|2007l ). in more detail, and for larger particle 
numbeJ^ 

What, otherwise, have we learnt from this study? The 
following points seem pertinent to us: 

• Analyses of violent relaxation usually neglect the ef- 
fects we have focussed on here, and assume that energy and 
mass of the final virialized system is equal to that of the 
initial configuration. This is true, for example, of the use 
of the s pherical collapse model in the context of cosmology 
(see e.g. ICoorav fc ShethI (|2002l )'). While it may be, in this 
context, a good approximation to neglect ejected mass — 
generically we would not expect to have the very large col- 
lapse factors observed here which are closely linked to this 
ejection — the non-trivial slow A^ dependences we have seen 
would urge some caution in this regard. Further in theoret- 
ical attempts to understand the properties of the virialized 
states, the same approximation is gener ically made. An evi - 
dent example is the famous proposal bv iLvnden-Belll (|l967t ) 
that the result of violent relaxation should be a state which 
maximizes a coarse-grained entropy in phase space subject 
to the constraints of the VP equations, and notably total 
energy and mass conservation. What we observe in our sim- 
ulations is that the density profile of the final state (and also 
velocity distributions, which we have not reported here) is 
very robust while the conservation of the energy and mass 
of the initial state appears to be of little relevance. This ob- 
servation may be important in attempts to understand the 
properties of these states, and in particular their dependence 
(or non-dependence) on initial conditions. 

• Our study shows that the assumption that the evolu- 
tion of such systems on these timescales is fully represented 
by the VP (or coUisionless) limit should also be treated with 
some caution. In our simulations we have found small, but 
significant, non-VP effects with a simple mass segregation 
test. In any case it is important to define clearly the appro- 



■^^ As the cost of numerical integration depends, essentially, on 

2° iBoilv et al. I 1I2OO2I ') gives results for a spherical collapse with the violence of the collapse, we reach our limit for the exact lattice 

10^ particles (but does not report the ejected mass and energy) at a particle number of order that reported (N=4143). 



Energy ejection in cold spherical collapse 19 



priate extrapolation to the VP limit, which we have seen 
may be less trivial than one might naively appears. We be- 
lieve the extrapolation we have defined here — in which 
particle number increases while the fluctuations above some 
scale are held fixed — is the appropriate one for this system. 

• Our results suggest that open self-gravitating systems 
may be able, theoretically, to eject an arbitarily large 
amount of energy on a dynamical time scale. This effect is 
reminiscent of the "gravothermal catastrophe", and indeed 
its origin is the same: because the gravitational potential is 
unbound below, an arbitrarily large kinetic energy can be 
gained by making some mass arbitrarily bound. The mech- 
anism, however, and timescale, are completely different. It 
is clear that the ejection we have seen is related to the col- 
lapse, which in turn is related to the spherical symmetry of 
the problem we have considered. What about more general 
initial conditions? We note that, although we expect smaller 
collapse factors for a generic initial condition, very large col- 
lapse factors may occur for non-symmetric initia l conditions. 
This q uestion has been considered at length in iBoilv et al.l 
l|2002t) . which presents a detailed numerical study of the gen- 
eralisation of the fluid SCM to axisymmetric and triaxial 
configurations. In the former singularities remain intact, and 
a relation Rmin oc A'^"^'^^ is found empirically to replace the 
Rmin OC A*'"^'^^ behavlour of the spherical case. In the triax- 
ial case the collapse factors are found to be typically finite, 
but they can be very large and no upper bound is placed on 
them. 

The study we have presented is a purely theoretical one 
of an idealized problem. It is interesting, of course, to con- 
sider whether, in particular, the "explosive" ejection of en- 
ergy we have found could have any direct relevance to real 
physical systems or realistic models of them in astrophysics 
or cosmology. For the former context it is relevant to con- 
sider the validity of the Newtonian limit which we have 
treated. To do so it is useful to write the ejected kinetic 
energy as 



(32) 



where and Mi are, respectively, the initial mass and 
(Newtonian) gravitational potential, and c is the speed of 
light. The first factor is, in fact, approximately the maximal 
value reached by the potential, and if this is small com- 
pared to unity we expect the Newtonian approximation to 
be valid. Thus our determination of would be expected 
to remain valid until this energy reaches of order the evident 
upper bound imposed by the rest energy of the ejected mass 
(which is of order that of the initial mass). We note also 
that if the initial radius of the sphere is Ri the Newtonian 
approximation can remain valid at all times provided the 
number of Poisson distributed masses satisfies 



f Ric 



\GMi 



(33) 



Thus, for example, if Mi and Ri are taken to be a solar 
mass and radius, respectively, the Newtonian approximation 
should be valid at all times during the collapse if A'' < 10^'', 
i.e., if the mass of the initially Poisson distributed "parti- 
cles" is greater than about 10~^* solar masses. Normalizing 
the ejected energy using our simulations above, one would 
obtain an ejected energy of about 10*® ergs for the largest 



value of A'^ we reported (A'' « 2 x 10^), and, extrapolating 
our results, an energy as large as that of a typical supernovae 
(~ lO^'^ ergs) for A'' ~ 10^^. Whether such initial conditions 
(of very uniformly distributed dark matter "particles" ) could 
possibly be produced in an astrophysical context and thus 
give rise to such purely gravitational "explosions", with ob- 
servational traces, is beyond the scope of our study. In the 
context of cosmology we note that the results we have found 
may be relevant, for example, in theories of structure forma- 
tion in "hot" or "warm" dark matter cosmologies. In these 
cases the initial spectrum of density fluctuations is cut-off 
abruptly below some cosmological scale, and the first struc- 
tures should form at such very large scales, without prior 
formation of structures at smaller scales. In this context the 
physical effects we have observed might then be expected to 
come into play, leading potentially also to implications for 
the numerical resolution provided by an A''-body discretisa- 
tion. 

Finally let us remark that rather than considering dif- 
ferent non-uniform and/or non-spherical initial conditions, 
the next natural step in this study is to consider the effects 
of the presence of initial velocity dispersion. While we would 
expect such dispersion to regulate the collapse and bound 
above the ejected mass, it would be interesting to see the 
full dependence of these quantities on A'^ and the amplitude 
of this velocity dispersion. It would be interesting to focus 
on a simple class of initial conditions for the velocities, e.g., 
"waterbag" distributions, studying very carefully the con- 
vergence of results for the virialized state as a function of 
A'^, and to characterize carefully also the velocity distribu- 
tions in the virialized states (which we have not discussed 
here) . Wi th respect to predict ions of theoretical models such 
as that of lLvnden-Belll l|l967^ . it would be interesting to ex- 
plore whether by appropriately modifying the energy /mass 
constraints on the final state better agreement may be ob- 
tained th an has been observed without such con s iderat ions 
(see, e.g., lArad fc JohanssonI (|2005h : lle^in et al.l (|2008h ). 
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and thus 5[R) = f{R)5o where 

fiR) = "27^ [v^(R - 3) - 3VT^ asin(%/r^)] . 

(A3) 

The behaviour of Eq. (|A3|I as _R 0, i.e., close to the col- 
lapse time, is 

f{R) = ^R~'^^ + 0{R-'/^) . (A4) 



APPENDIX A: GROWTH OF LINEARIZED 
PERTURBATIONS 

After some simple algebra, Eq. (O reduces to 

^dPS f I \ dS f 3 4\ 3 , „ ,,,, 
'dR^[R-') + dR[R^~R)-R^'-'^ 

where we have taken R{t = 0) = 1. The general solution can 
be written as 

^2-^ (yR{R -3)+ SVR^ a8mh{^/R~T)J 
Using the appropriate initial conditions for our case, 

it follows that 

Ai=0 ,A2 = -So/2, 



